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ABSTRACT 

Many  problems  in  physics  and  engineering  can  be  reduced  to  the  Sturm- 
Liouville  (  S-L )  Problem.  For  some  classes  of  the  S-L  problem  we  develop 
•  a  method  that  can  rapidly  result  in  computational  solutions.  Let  us  assume 

that  we  can  write  the  problem  as  follows: 

d2U.(r) 

-^-  +  (K(r)-X,)Uj(r)  =  0 

where  L  is  an  eigenvalue  and  K(r)  is  some  well  behaved  (i.e.  positive 
definite)  function  of  r.  Now  let  k  correspond  to  some  average  value  of  K(r) 
,  over  the  domain  of  K(r).  The  boundary  conditions  are  general.  Then  we 

easily  solve  the  following  problem: 

d2Vj(r) 

“l^  +  [k-Wr)  =  0. 

This  yields  global  basis  functions.  We  determine  an  efficient  method 
'  to  obtain  solutions  of  the  form: 
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in  which  we  determine  the  best  expansion  coefficients  a-  as  well  as  the 
eigenvalues  Lj’s  in  terms  of  the  Ij’s  and  the  a^’s.  We  also  develop  several 
perturbation  methods  from  the  technique.  The  method  is  then  applied  to 
typical  problems  found  in  acoustics  and  quantum  mechanics. 

I.  INTRODUCTION 

One  of  the  most  interesting  problems  of  the  Sturm-Liouville  class  pertains 
to  the  solution  of  the  vertical  component  of  the  normal  mode  of  sound  in 
a  wave  guide.1, 2  We  shall  discuss  a  general  method,  based  on  separation  of 
variables,  which  is  ideally  suited  for  a  range-independent  wave  guide  but 
has  been  extended  to  certain  range-dependent  problems  by  mode  coupling.3, 4 
It  is  safe  to  say  that  this  appealing  method  has  been  explored  by  a  large 
number  of  researchers,  and  several  numerical  strategies  have  been  advanced 
to  find  solutions.  Generally,  one  is  required  to  solve  the  vertical  component 
of  the  problem  subject  to  a  variety  of  environmental  constraints.  The  range 
part  of  the  solution  is  analytical.  The  vertical  solutions  are,  in  fact,  eigen¬ 
functions  that  fall  into  the  category  of  a  Sturm-Liouville  problem.5  For  the 
isovelocity  case  the  eigenvalues  for  a  layered  bottom  and  a  pressure  release 
surface  are  simple  to  find,  and  the  solutions  in  the  water  column  are  simple 
sine  function.  For  variable  velocity  profiles  the  water  column  can  be 
represented  as  a  sum  of  isovelocity  layers  that  have  discrete  changes  at  the 
isovelocity  interface  or  the  problem  can  be  solved  numerically.  For  certain 
classes  of  velocity  profiles  (with  linear  behavior)  the  Airy  function,  which 
seems  likely  to  be  one  of  the  fastest  of  the  techniques,  can  be  employed. 

•tr 

One  may  ask  the  question,  “So  why  another  method?”  The  principal 
reason  is,  that  for  certain  applications  it  is  desirable  to  obtain  the  normal 
mode  solution  in  closed  analytical  form.  The  motivation  behind  this  research 
is  to  formulate  the  normal  mode  solution  in  such  a  manner  that  the  solution 
can  be  written  in  the  form  of  a  spherical  representation.  This  formulation 
is  required  for  some  problems  pertaining  to  scattering  from  objects  in  a 
waveguide.  It  is  quite  easy  to  see  that  the  isovelocity  case  is  just  such  a 
representation.  This  fact  suggests  that  a  perturbation  about  the  isovelocity 
solution  be  attempted.  However,  conventional  perturbation  theory  is  too 
limiting  r  d  works  only  with  very  small  departures  from  the  isovelocity 
case. 

The  purpose  of  this  work  is  to  develop  a  method  using  the  isovelocity 
solution  as  a  set  of  basic  functions  that  span  the  solution  space  that  adequately 
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represents  the  exact  eigenvalue  problem.  This  method  uses  Sturm-Liouville 
theory  and  completeness.  A  nice  outgrowth  of  the  method  to  be  developed 
is  to  show  why  the  old  perturbation  theory  is  limited  and  how  to  extend  it 
by  adding  a  simple  term,  so  that  the  new  perturbation  method  will  work  for 
more  severe  variations  from  the  isovelocity  case. 


II.  THEORETICAL  DEVELOPMENTS 


An  updated  perturbation  method 

We  write  the  exact  solution  for  unperturbed  cases  as  follows: 


d2Vi 

dz2 


+ 


ko  ~  Vi  -  0 


and  the  desired  solution  as  follows: 


d2Uj 

dz2 


k2(z)  -  X; 


Ui  =0. 


For  convenience,  rewrite  Equation  2  as  follows: 


d2U;  t 
dz2  + 


Uj  =  QUj  ( 


0) 


(2) 


(3) 


where  Q  =  1^  -  k2(z)  +  Xi  -  X?  =  q  +  A  Xj  and  where  q(Z)  =  k^  -  k2(z)  and 
AX;  =Xj-X|V 


For  the  isovelocity  case  we  can  rewrite  Equation  1  as 


dVi 

dz2 


+  af  Vi 


(4) 


where  af  »  k£  -  xf- 

Impose  the  following  orthonormality  conditions: 


=  5 


»i  and 


l 

P 


Ui,-Uj 


-8i 


(5) 
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Assume  closure  so  that  we  can  express: 


N 

Ui=  I  ajj  Vj . 

j  =  i 


(6) 


Insert  Equation  6  into  Equation  5  and  make  use  of  the  Galerkin  weighted 
residual  method  we  arrive  at 


N  /  2  2)  N 

X  aij  CX;  -  aj  Vj  =  I  ajj  (q(z)  +  A  XjJvj , 
i  =  l  \  I  j  =  l 

where  i  =  1,  2,  3,  . . .  N. 

By  integrating  the  overlap  of  Equation  7  with  V;  we  obtain: 

o  N 

X;  »  X;  —  q;;  -  X  ajj  q*j  > 

a-  i,tj 

-  *11 

where  aij  =  — \  This  yields  the  Eigenvalue  Correction  Equation. 
aii 

Integrate  the  overlap  of  Equation  7  with  yk  to  obtain: 


aik  |®i  -  ctjj  -  A  Xk]  -  X  Rjk  ajj  -  qik  t 

'  '  j  *  i 


where  k  =  1 ,  2,  3,  . . .  N  and  i  *  k 


(7) 


(8) 


(9) 


This  is  an  eigenvalue  equation.  We  can  use  the  expansion  for  the 
eigenvalue  above  to  rewrite  the  equation  as  follows: 


iik  |af  -  al  +  Hk|  -  X  <fck  h  =  q*  , 

j*(i.k) 


(10) 


where  k  =  1,  2,  3  . . .  N,  i  *  k,  and  the  diagonal  terms  are. 

This  expression  will  prove  useful  later.  The  Hk’s  are  the  higher  order 
terms  and  are  negligible  in  many  cases.  The  diagonal  terms  are  almost 
always  greater  than  the  off-diagonal  terms.  As  a  first  iteration  via  Gauss- 
Seidel  we  obtain 


&ik  = 


q.k 


q.k 


2  2  2  2 
<*j  -  ak  +  qkk  -  A  Xk  a;  -  ak  -  Hk 


(11) 
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where 

N 

Hi  =  ~  X  aij  Qij  t 
j*> 


(12) 


The  term  Hk  can  be  a  big  improvement  over  the  ordinary  perturbation 
term  in  that  often  Hk  >  q^.  Thus,  we  expect  that  generally  we  will  have 
convergence  if  the  master  matrix  (defined  by  Eq.  10)  is  diagonally  strong. 
Further,  this  expression  shows  that  the  old  theory  overestimated  the  expansion 
coefficient. 

The  new  complete  perturbation  expansion  for  the  eigenvalue  is  now 


i*j  <Xj  -  dj  +  Hj 


This  expression  indicates  that  the  first-order  correction  is  still  the  same  but 
that  the  higher-order  corrections  are  over/underestimated  in  the  old  theory. 

A  solution  to  the  eigenvalue  problem  and  an  eigenvalue  approach 

Reorder  the  above  equations  selectively  to  arrive  at 

Qa  =  la  (14) 

where  Q;j  =  -  qij  i  *  j  and  Q;;  =  X?  -  qu- 

We  thus  have  an  eigenvalue  problem  for  the  l’s.  In  Equation  14  the 
matrix  Q  is  real,  symmetric,  and  diagonally  strong.  A  rather  good  strategy 
to  solve  this  eigenvalue  problem  is  to  use  the  Householder  method  to 
reduce  Q  to  tridiagonal  form.  For  an  n  x  n  matrix,  this  only  takes  n-2 
orthogonal"  transformations.  Since  the  eigenvalues  of  this  very  stable  and 
fast  reduction  method  are  unchanged,  we  need  only  to  find  the  eigenvalues 
of  a  tridiagonal  matrix,  which  can  be  done  in  n  operations.  We  have  also 
tried  the  more  general  QR  algorithm  for  the  above  problem;  although  we 
arrive  at  the  same  set  of  eigenvalues,  the  method  is  understandably  somewhat 
slower.  The  quantity  a  in  the  above  expression  is  rotational  (unitary  for 
complex  eigenvalues)  and  should  correspond  to  the  expansion  coefficients 
required  in  reconstructing  the  exact  eigenvalues.  Indeed,  the  fact  that  matrix  a 
is  rotational  or  unitary  under  certain  conditions  guarantees  that  the  newly 
constructed  eigenvalues  are  orthogonal,  as  is  required  for  Sturm-Liouville 
problems. 

However,  there  are  some  difficulties  when  one  tries  to  retain  a  from 
the  process  using  either  the  QL  or  QR  algorithms.  One  problem  is  that  the 
relative  phases  between  the  eigenfunctions  are  not  retained,  which  leads  to 


-«V-;  wf ; 


sxy'S;. 
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erroneous  results  when  summing  the  normal  mode  series.  The  second  problem 
appears  to  be  the  accumulation  of  numerical  errors  for  large  numbers  of 
normal  modes.  Thus,  we  use  the  above  method  just  to  obtain  the  eigenvalues 
and  revert  to  the  master  equation  that  defines  a  along  with  Hk  to  obtain  the 
expansion  coefficients.  This  method  seems  to  work  well.  Interestingly,  the 
new  perturbation  method  yields  rather  close  agreement  with  the  exact  method 
for  fairly  strong  perturbations,  while  the  old  perturbation  method  performs 
poorly.  The  next  section  presents  examples. 

III.  NUMERICAL  EXAMPLES 

For  the  example  we  chose  a  pressure  release  surface  with  a  water 
column  of  300  m  over  a  fluid  half-space.  The  velocity  profile  (Fig.  1) 
varies  from  1510  m/s  at  the  surface  to  1476  m/s  at  100  m  linearly  to 
1520  m/s  at  the  bottom.  The  fluid  half-space  has  a  compressional  speed  of 
1600,  a  density  of  1.6,  and  an  attenuation  of  0.5  dB/m.  Absorption  in  the 
bottom  is  treated  by  the  standard  perturbation  treatment. 


Figure  1.  Diagram  of  the  waveguide  configuration  and 
parameters  used  in  the  example. 
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We  try  several  strategies  for  obtaining  the  isovelocity  spanning  eigen¬ 
functions.  The  wisest  and  safest  strategy  is  to  choose  the  isovelocity  to 
correspond  to  the  smallest  value  of  the  velocity  profile,  in  this  case,  1476. 
This  value  guarantees  that  we  will  not  miss  any  modes.  It  is  all  right  to 
span  the  solution  space  with  a  larger  number  of  isovelocity  eigenfunctions. 
The  boundary  condition  at  the  water-ocean  bottom  interface  will  exclude 
imaginary  modes  if  they  are  encountered.  One  might  also  choose  the 
isovelocity  as  the  average  value  along  the  water  column.  The  advantage 
here  would  be  to  reduce  errors  in  the  predicted  eigenvalues,  but  such  a 
choice  obviously  might  exclude  the  highest  order  modes.  As  a  final  strategy 
we  choose  a  small  isovelocity  value  of  1400  m/s  to  determine  the  robust¬ 
ness  of  the  method  and  to  see  whether  or  not  we  could  account  for  some 
of  the  continuum  which  we  have  ignored.  This  last  strategy  did  not  yield 
a  correct  answer. 

We  also  compared  our  results  with  Michael  D.  Collins’  finite-element 
normal  mode  code.  We  employed  a  Taylor’s  expansion  of  the  inverse  of 
the  variable  velocity  about  the  isovelocity  and  dropped  all  but  the  linear 
terms  to  obtain  analytical  integrals.  This  approach  introduced  an  error  of 
a  few  parts  per  thousand,  so  that  a  comparison  with  other  methods  should 
only  be  accurate  to  about  three  places.  Table  1  compares  the  seven  eigen¬ 
values  for  a  50-Hz  signal  for  the  low  and  the  average  reference  velocities 
to  those  predicted  by  Collins’  code.  The  two  reference  velocities  yield 
essentially  exact  agreement.  Agreement  with  the  Collins  value  is  within 
the  accuracy  tolerance  of  a  few  pans  per  thousand.  We  also  list  the  eigenvalues 
for  the  isovelocity  case  using  the  lower  reference  number.  A  comparison 
of  the  eigenfunction  using  the  two  reference  wave  numbers  are  in  very 
good  agreement.  Figure  2  illustrates  the  seven  modes  constructed  using  the 
lower  wave  number.  It  is  clear,  as  expected,  that  channeling  occurs  at 
100  m.  We  ran  the  code  for  this  problem  to  400  Hz,  and  58  modes  and  the 
results  were  apparently  successful. 


Table  1 .  Eigenvalues  for  waveguide  problem  at  50  Hz. 


Mode 

FEMOOE 

(M.Collins) 

NMEIG 

(Lower) 

NMEIG 

(Ave) 

NMEIG 

(ISO) 

1 

0.2113 

0.2116 

0.2116 

0.2126 

0.2091 

0.2086 

0.2087 

0.2119 

0.2072 

0.2069 

0.2070 

0.2107 

0.2054 

0.2054 

0.2055 

0.2092 

5 

0.2034 

0.2032 

0.2033 

0.2070 

6 

0.2008 

0.2007 

0.2006 

0.2044 

7 

0.1978 

0.1974 

0.1974 

0.2012 
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We  now  want  to  compare  the  exact  solution  with  both  the  old  and  new 
perturbation  theories.  In  all  cases  we  use  the  lower  reference  isovelocity  to 
generate  results.  A  good  test  of  the  method  is  to  examine  the  expansion 
coefficients  for  the  three  cases  (Tables  2,  3,  and  4).  The  exact  expansion 
coefficients  (Table  2)  and  those  from  the  old  perturbation  theory  (Table  3) 
are  widely  different,  but  the  coefficients  from  the  new  method  agree  reasonably 
well  with  the  exact  method.  A  better  idea  of  the  agreement  can  be  seen  by 
comparing  the  exact  eigenvalues  (Fig.  2)  with  those  produced  by  the  new 
method  illustrated  in  Figure  3.  Visual  agreement  is  also  quite  good.  The 
old  perturbation  theory  does  not  even  reproduce  the  correct  modal  ordering. 


Table  2.  Normal  mode  expansion  coefficients  (exact). 


Mode 

1 

2 

3 

4 

5 

6 

7 

8 

1 

1.000 

0.384 

-0.102 

-0.241 

-0.151 

0.014 

0.086 

0.056 

2 

0.258 

1.000 

0.126 

-0.168 

-0.126 

-0.035 

-0.031 

0.039 

3 

-0.095 

0.129 

1.000 

0.133 

-0.134 

-0.095 

-0.028 

0.006 

4 

-0.121 

-0.162 

0.128 

1.000 

0.149 

-0.124 

-0.106 

-0.032 

5 

-0.680 

-0.116 

-0.127 

0.148 

1.000 

0.129 

-0.134 

-0.119 

6 

0.001 

-0.036 

-0.094 

0.128 

0.134 

1.000 

0.111 

-0.110 

7 

0.034 

0.029 

-0.028 

-0.111 

-0.141 

-0.111 

1.000 

0.119 

Table  3.  Normal  mode  expansion  coefficients  (old  perturbation). 


Mode 

1 

2 

3 

4 

5 

6 

7 

8 

1 

1.000 

2.085 

-0.402 

-0.102 

-0.050 

0.0005 

0.005 

0.0006 

2 

-2.085 

1.000 

0.928 

-0.367 

-0.072 

-0.029 

0.002 

0.008 

3 

0.402 

-0.928 

1.000 

0.652 

-0.259 

-0.053 

-0.018 

-0.007 

4 

0.102 

0.367 

-0.652 

1.000 

0.510 

-0.198 

-0.054 

-0.010 

5 

0.050 

0.072 

0.259 

-0.510 

1.000 

0.383 

-0.158 

-0.062 

6 

-0.0005 

0.029 

0.053 

0.198 

-0.383 

1.000 

0.287 

-0.114 

7 

-0.005 

-0.002 

0.018 

0.054 

0.158 

-0.287 

1.000 

0.217 

Table  4.  Normal  mode  expansion  coefficients  (new  perturbation). 


Mode 

1 

2 

3 

4 

5 

6 

7 

8 

1 

1.000 

0.338 

-0.164 

-0.080 

-0.063 

0.0009 

0.013 

0.002 

2 

0.258 

1.000 

0.163 

-0.157 

-0.055 

-0.034 

0.003 

0.019 

3 

-0.124 

0.160 

1.000 

0.154 

-0.142 

-0.049 

-0.052 

-0.013 

4 

-0.056 

-0.151 

0.151 

1.000 

0.155 

-0.135 

-0.062 

-0.017 

5 

-0.046 

-0.052 

-0.138 

0.154 

1.000 

0.143 

-0.131 

-0.085 

6 

0.0006 

-0.033 

-0.048 

-0.135 

0.145 

1.000 

0.131 

-0.114 

7 

0.010 

0.003 

-0.025 

-0.061 

-0.131 

0.129 

1.00 

0.117 
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As  frequency  increases,  channeling  becomes  more  pronounced  for  the 
lower  modes.  This  seems  to  create  some  problem  with  the  new  perturbation 
method  because  at  that  point  the  governing  matrix  begins  to  depart  from 
a  diagonally  strong  matrix  so  that  Gauss-Seidel  cannot  be  used  in  a  first 
iteration.  However,  a  second  or  third  iteration  may  be  adequate  for  that 
case,  and  it  could  be  done  rapidly  by  computer.  This  problem  existed  for 
only  the  two  lowest  modes  at  400  Hz.  The  higher  modes  agreed  quite  well 
for  the  new  perturbation  method  and  the  exact  method. 


IV .  COMMENTS 

The  expansion  method  presented  here  is  useful  for  formulating  a  normal- 
mode  solution  in  a  closed  mathematical  form.  An  outgrowth  of  the 
development  was  a  new  perturbation  theory  that  is  more  powerful  than  an 
earlier  one  commonly  used  in  quantum  mechanics.  The  method  is  quite 
fast,  however,  and  calculated  58  modes  in  about  two  minutes  on  a  Vax  8650. 
With  further  development  and  testing,  it  may  prove  to  be  a  very  useful 
general  method.  Another  advantage  of  the  method  would  be  its  use  in  a 
range-dependent  normal-mode  development. 

It  is  also  possible  to  develop  a  deep-water  model  by  matching  the 
expansion  wave  functions  at  some  point  in  the  water  column  (when  the 
velocity  profile  has  a  constant  slope)  to  an  Airy  function. 
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